1 Giotto Suite Object - A Flexible Spatial Analysis Framework

The Giotto object holds spatial-omic data and facilitates its analysis and visualization. Central to the design of the Giotto Suite object is the concept of spatial units and feature types. Most spatial data used to be pre-aggregated, with the expression information ‘belonging’ to an immutable unit of study. This lent themselves to similar analysis and organization as single-cell-based analyses. With recent methods, raw feature data is now decoupled from the morphological information, opening many doors for how the data can be analyzed and combined.

A spatial unit (spat_unit) is a set of polygonal information that can be used to define what you want your functional unit to be when performing analysis. This can be at different spatial sizes such as polygons that only select the nuclear region, the entire cell, or even larger tissue structures.

A feature type (feat_type) is which modality of feature you are analyzing such as rna or protein.

Analyses should then be organized first by which spat_unit is being studied, and then by which feat_type is being used to understand that spat_unit. Giotto’s updated object adopts a nested organization that reflects this in order to maximize flexiblity and ease of exploration.

2 Setup

i_p = installed.packages()
if(!"Giotto" %in% i_p)
  devtools::install_github("drieslab/Giotto@suite")
# Ensure Giotto Data is installed
if(!"GiottoData" %in% i_p)
  devtools::install_github("drieslab/GiottoData")


library(data.table)
library(Giotto)
library(GiottoData)

# Ensure the Python environment for Giotto has been installed
genv_exists = checkGiottoEnvironment()
if(!genv_exists){
  # The following command need only be run once to install the Giotto environment
  installGiottoEnvironment()
}

3 Creating a Giotto Object

Assembling a complete object for spatial analysis requires expression and spatial information.
This can be provided in the form of aggregated information (expression matrices with paired spatial locations) or as raw feature data (feature detections or images) and polygonal information.

This tutorial describes how you can assemble a Giotto object starting from subcellular data using a subset of the Vizgen MERSCOPE Mouse Brain Receptor Map data release.

3.1 Load Data

## provide path to vizgen folder
data_path = system.file('/Mini_datasets/Vizgen/Raw/',
                        package = 'GiottoData')

## 0.1 path to transcripts ####
# --------------------------- #
## each transcript has x, y and z coordinate
tx_path = paste0(data_path, '/', 'vizgen_transcripts.gz')
tx_dt = data.table::fread(tx_path)
tx_dt = tx_dt[,.(global_x, -global_y, gene, global_z)]

## 0.2 path to cell boundaries folder ####
# -------------------------------------- #
## vizgen already provides segmentation information in .hdf5 files
## Here I have already converted the hdf5 files to a simple data.table format

boundary_path = paste0(data_path, '/cell_boundaries/')

z0_polygon_DT = fread(paste0(boundary_path, '/', 'z0_polygons.gz'))
z1_polygon_DT = fread(paste0(boundary_path, '/', 'z1_polygons.gz'))

## 0.3 path to images ####
# ---------------------- #
DAPI_z0_image_path = paste0(data_path, '/', 'images/mini_dataset_dapi_z0.jpg')
DAPI_z1_image_path = paste0(data_path, '/', 'images/mini_dataset_dapi_z1.jpg')

3.2 New Object Types

3.2.1 giottoPolygon

# giottoPolygon creation from data.table
z0_polygons = createGiottoPolygonsFromDfr(name = 'z0',
                                          segmdfr = z0_polygon_DT)
z1_polygons = createGiottoPolygonsFromDfr(name = 'z1',
                                          segmdfr = z1_polygon_DT)

z0_polygons
## An object of class giottoPolygon
## spat_unit : "z0"
## Spatial Information:
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 498, 1  (geometries, attributes)
##  extent      : 6399.244, 6903.243, -5152.39, -4694.868  (xmin, xmax, ymin, ymax)
##  coord. ref. :  
##  names       :                                 poly_ID
##  type        :                                   <chr>
##  values      :  40951783403982682273285375368232495429
##                240649020551054330404932383065726870513
##                274176126496863898679934791272921588227
##  centroids   : NULL
##  overlaps    : NULL
plot(z0_polygons)

3.2.2 giottoPoints

# giottoPoints creation from data.table
tx = createGiottoPoints(tx_dt, feat_type = 'rna')

tx
## An object of class giottoPoints
## feat_type : "rna"
## Feature Information:
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 80343, 3  (geometries, attributes)
##  extent      : 6400.037, 6900.032, -5149.983, -4699.979  (xmin, xmax, ymin, ymax)
##  coord. ref. :  
##  names       : feat_ID global_z feat_ID_uniq
##  type        :   <chr>    <int>        <int>
##  values      :    Mlc1        0            1
##                 Gprc5b        0            2
##                   Gfap        0            3
plot(tx)

3.2.3 giottoLargeImage

ultra_mini_extent = ext(c(6400.029, 6900.037, -5150.007, -4699.967 ))
image_paths = c(DAPI_z0_image_path, DAPI_z1_image_path)
image_names = c('dapi_z0', 'dapi_z1')

imagelist = createGiottoLargeImageList(raster_objects = image_paths,
                                       names = image_names,
                                       extent = ultra_mini_extent)

imagelist[[1]]
## An object of class giottoLargeImage : "dapi_z0"
## Image extent            : 6400.029, 6900.037, -5150.007, -4699.967 (xmin, xmax, ymin, ymax)
## Original image extent   : 6400.029, 6900.037, -5150.007, -4699.967 (xmin, xmax, ymin, ymax)
## Scale factor            : 0.108626547903541, 0.10862659908279 (x, y)
## Resolution              : 9.2058527063567, 9.20584836903386 (x, y)
## Layers                  : 1 
## Estimated max intensity : 208 
## Estimated min intensity : 0 
## Values                  : integers
## File path               : '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/GiottoData//Mini_datasets/Vizgen/Raw//images/mini_dataset_dapi_z0.jpg'
plot(imagelist[[1]])

3.3 Object Creation

Giotto objects are created through a call to createGiottoObject(). Inputs to to this function should be supplied as named lists of data where the names are used to determine the data nesting they should be assigned as.

viz = createGiottoObject(feat_info = list('rna' = tx_dt),
                         spatial_info = list('z0' = z0_polygons))
instructions(viz, 'return_plot') = FALSE # set instructions to prevent returning of plots

4 Subcellular Workflow

The next steps walk through the usual method of processing a Giotto object that has been assembled using subcellular information

4.1 Add Spatial Locations

Generate centroids to use as convenient spatial representations of where your cells are. This information gets stored in the @spatial_locs slot

viz = addSpatialCentroidLocations(gobject = viz,
                                  poly_info = 'z0',
                                  return_gobject = TRUE)
# visualize
spatPlot2D(viz)

4.2 Add Images

viz = addGiottoImage(viz, largeImages = imagelist[1])
## Warning in spatInSituPlotPoints(viz, show_image = TRUE, largeImage_name = "dapi_z0", : You need to select features (feats) and modify feature types (feat_type) if you want to show individual features (e.g. transcripts)
## select image done
## plot image layer done
## plot polygon layer done

4.3 Aggregate Expression Matrix

Using the polygonal information, sum up all of the transcript point detections that fall within. This is can then be used as an expression count matrix. This information is then stored in the @expression slot.

# Calculate transcripts in 'rna' overlapped
# by the specified polygons in 'z0' (aggregation)
viz = calculateOverlapRaster(viz,
                             spatial_info = 'z0',
                             feat_info = 'rna')

# Send overlaps information to the expression slot as a new expression matrix
# of spat_unit 'z0' and feat_type 'rna' with a name of 'raw'
viz = overlapToMatrix(viz,
                      poly_info = 'z0',
                      feat_info = 'rna',
                      name = 'raw')
## Warning in spatInSituPlotPoints(viz, show_polygon = TRUE, polygon_feat_type = "z0", : You need to select features (feats) and modify feature types (feat_type) if you want to show individual features (e.g. transcripts)
## plot polygon layer done

From here on, the normal aggregate type analyses can continue.

4.4 Exploring the Object

We have added several objects to the Giotto object now. Directly returning the object gives a summary of what information is contained within.
For a more detailed look at a specific slot and type of data, Giotto has several show functions.

# Returning the object provides a summary
viz
## An object of class giotto 
## >Active spat_unit:  z0 
## >Active feat_type:  rna 
## [SUBCELLULAR INFO]
## polygons      : z0 
## features      : rna 
## [AGGREGATE INFO]
## expression -----------------------
##   [z0][rna] raw
## spatial locations ----------------
##   [z0] raw
## attached images ------------------
## giottoLargeImage : dapi_z0 
## 
## 
## Use objHistory() to see steps and params used
# Closer look at @expression slot contents
showGiottoExpression(viz)
## └──Spatial unit "z0"
##    └──Feature type "rna"
##       └──Expression data "raw" values:
##             An object of class exprObj : "raw"
##             spat_unit : "z0"
##             feat_type : "rna"
##             provenance: z0 
##             
##             contains:
##             559 x 498 sparse Matrix of class "dgCMatrix"
##                                                     
##             Mlc1   . . . . . . . .  1 . . 1 . ......
##             Gprc5b . . 1 . 3 . . .  2 . 4 2 . ......
##             Gfap   . . . 2 1 1 . . 48 . 1 1 . ......
##             
##              ........suppressing 485 columns and 553 rows 
##                                                      
##             Chat     . . . . . . . . . . . . . ......
##             Blank-58 . . . . . . . . . . . . . ......
##             Xcr1     . . . . . . . . . . . . . ......
##             
##              First four colnames:
##              40951783403982682273285375368232495429
##              240649020551054330404932383065726870513
##              274176126496863898679934791272921588227
##              323754550002953984063006506310071917306 
## 

5 Multiple Z-Stacks and Spatial Unit Provenance

What about when the data is presented as multiple confocal z-stacks? The Vizgen MERSCOPE Mouse Brain Receptor Map dataset is actually one such example, with 7 layers of different polygons alongside matching transcript detections at each Z layer. This data should not be treated as fully 3D since the z-layers are close enough together that they sample multiple times from each individual cell as opposed to entirely different cells in a different z layer.
This problem illustrates another way in which spat_unit can be used in Giotto. Giotto handles this by defining each z-layer as its own spat_unit and then performing separate aggregations at each z-layer. These spat_units are then merged them into a single combined spat_unit for analysis.
Since the resulting combined spat_unit is different depending on which layers were used, that information is tracked as provenance.

5.1 Remake Giotto Object

In light of this, let’s remake our object.

# create with 2 spatial units
viz = createGiottoObject(feat_info = list('rna' = tx),
                         spatial_info = list('z0' = z0_polygons,
                                             'z1' = z1_polygons))
instructions(viz, 'return_plot') = FALSE # set instructions to prevent returning of plots

# generate spatial locations
viz = addSpatialCentroidLocations(gobject = viz,
                                  poly_info = paste0('z',0:1),
                                  provenance = list('z0', 'z1'),
                                  return_gobject = TRUE)

5.2 Redo Aggregation

5.2.1 Layer-Specific Aggregations

viz = calculateOverlapRaster(viz,
                             spatial_info = 'z0',
                             feat_info = 'rna',
                             feat_subset_column = 'global_z',
                             feat_subset_ids = 0)

viz = overlapToMatrix(viz,
                      poly_info = 'z0',
                      feat_info = 'rna',
                      name = 'raw')

viz = calculateOverlapRaster(viz,
                             spatial_info = 'z1',
                             feat_info = 'rna',
                             feat_subset_column = 'global_z',
                             feat_subset_ids = 1)

viz = overlapToMatrix(viz,
                      poly_info = 'z1',
                      feat_info = 'rna',
                      name = 'raw')

5.2.2 Collapse Into Common ‘aggregate’ spat_unit

# combine the calculated data for z1 and z0 into a new spatial unit called
# 'aggregate'
viz = aggregateStacks(gobject = viz,
                      spat_units = c('z0', 'z1'),
                      feat_type = 'rna',
                      values = 'raw',
                      summarize_expression = 'sum',
                      summarize_locations = 'mean',
                      new_spat_unit = 'aggregate')
## Warning in spatInSituPlotPoints(viz, show_polygon = TRUE, polygon_feat_type = "aggregate", : You need to select features (feats) and modify feature types (feat_type) if you want to show individual features (e.g. transcripts)
## plot polygon layer done

From here the usual aggregate analyses can be performed after filtering and normalization

viz = filterGiotto(viz,
                   expression_threshold = 1,
                   feat_det_in_min_cells = 1,
                   min_det_feats_per_cell = 1)
viz = normalizeGiotto(viz)
## 
## first scale feats and then cells
# You can see that there are now z0, z1, and 'aggregate' as sets of spatial units
viz
## An object of class giotto 
## >Active spat_unit:  aggregate 
## >Active feat_type:  rna 
## [SUBCELLULAR INFO]
## polygons      : z0 z1 aggregate 
## features      : rna 
## [AGGREGATE INFO]
## expression -----------------------
##   [z0][rna] raw
##   [z1][rna] raw
##   [aggregate][rna] raw normalized scaled
## spatial locations ----------------
##   [z0] raw
##   [z1] raw
##   [aggregate] raw
## 
## 
## Use objHistory() to see steps and params used
showGiottoExpression(viz)
## ├──Spatial unit "z0"
## │  └──Feature type "rna"
## │     └──Expression data "raw" values:
## │           An object of class exprObj : "raw"
## │           spat_unit : "z0"
## │           feat_type : "rna"
## │           provenance: z0 
## │           
## │           contains:
## │           494 x 498 sparse Matrix of class "dgCMatrix"
## │                                                   
## │           Mlc1   . . . . . . . .  . . . . . ......
## │           Gprc5b . . 1 . 1 . . .  1 . 2 . . ......
## │           Gfap   . . . 1 1 . . . 27 . . . . ......
## │           
## │            ........suppressing 485 columns and 488 rows 
## │                                                    
## │           Qrfpr    . . . . . . . . . . . . . ......
## │           Chat     . . . . . . . . . . . . . ......
## │           Blank-58 . . . . . . . . . . . . . ......
## │           
## │            First four colnames:
## │            40951783403982682273285375368232495429
## │            240649020551054330404932383065726870513
## │            274176126496863898679934791272921588227
## │            323754550002953984063006506310071917306 
## │        
## ├──Spatial unit "z1"
## │  └──Feature type "rna"
## │     └──Expression data "raw" values:
## │           An object of class exprObj : "raw"
## │           spat_unit : "z1"
## │           feat_type : "rna"
## │           provenance: z1 
## │           
## │           contains:
## │           494 x 504 sparse Matrix of class "dgCMatrix"
## │                                                   
## │           Mlc1   . . . . . . . . . . .  1 . ......
## │           Gprc5b . . . . . . . 2 . . .  1 . ......
## │           Gfap   . 3 . 2 . 1 4 . . . . 18 . ......
## │           
## │            ........suppressing 491 columns and 488 rows 
## │                                                    
## │           Qrfpr    . . . . . . . . . . . . . ......
## │           Chat     . . . . . . . . . . . . . ......
## │           Blank-58 . . . . . . . . . . . . . ......
## │           
## │            First four colnames:
## │            40951783403982682273285375368232495429
## │            17685062374745280598492217386845129350
## │            223553142498364321238189328942498473503
## │            240649020551054330404932383065726870513 
## │        
## └──Spatial unit "aggregate"
##    └──Feature type "rna"
##       ├──Expression data "raw" values:
##       │     An object of class exprObj : "raw"
##       │     spat_unit : "aggregate"
##       │     feat_type : "rna"
##       │     provenance: z0 z1 
##       │     
##       │     contains:
##       │     494 x 478 sparse Matrix of class "dgCMatrix"
##       │                                             
##       │     Mlc1   . . . . . .  1 . . 1 . . . ......
##       │     Gprc5b . 1 . 3 . .  2 . 4 2 . 1 . ......
##       │     Gfap   2 . 2 1 . . 45 . 1 1 . . . ......
##       │     
##       │      ........suppressing 465 columns and 488 rows 
##       │                                              
##       │     Qrfpr    . . . . . . . . . . . . . ......
##       │     Chat     . . . . . . . . . . . . . ......
##       │     Blank-58 . . . . . . . . . . . . . ......
##       │     
##       │      First four colnames:
##       │      240649020551054330404932383065726870513
##       │      274176126496863898679934791272921588227
##       │      323754550002953984063006506310071917306
##       │      87260224659312905497866017323180367450 
##       │  
##       ├──Expression data "normalized" values:
##       │     An object of class exprObj : "normalized"
##       │     spat_unit : "aggregate"
##       │     feat_type : "rna"
##       │     provenance: z0 z1 
##       │     
##       │     contains:
##       │     494 x 478 sparse Matrix of class "dgCMatrix"
##       │                                                                                             
##       │     Mlc1    .       .        .       .        . .  5.537748 . .        6.080373 . .        .
##       │     Gprc5b  .       6.658211 .       7.699722 . .  6.522136 . 7.927329 7.069674 . 5.338912 .
##       │     Gfap   10.74423 .        8.16347 6.128572 . . 10.998911 . 5.945000 6.080373 . .        .
##       │                  
##       │     Mlc1   ......
##       │     Gprc5b ......
##       │     Gfap   ......
##       │     
##       │      ........suppressing 465 columns and 488 rows 
##       │                                              
##       │     Qrfpr    . . . . . . . . . . . . . ......
##       │     Chat     . . . . . . . . . . . . . ......
##       │     Blank-58 . . . . . . . . . . . . . ......
##       │     
##       │      First four colnames:
##       │      240649020551054330404932383065726870513
##       │      274176126496863898679934791272921588227
##       │      323754550002953984063006506310071917306
##       │      87260224659312905497866017323180367450 
##       │  
##       └──Expression data "scaled" values:
##             An object of class exprObj 
##             for spatial unit: "aggregate" and feature type: "rna"
##               Provenance: z0 z1
##             
##             contains:
##             494 x 478 dense matrix of class "dgeMatrix"
##             
##                          [,1]       [,2]       [,3]       [,4]
##             Mlc1   -0.7655816 -0.5601177 -0.4625350 -0.4286685
##             Gprc5b -1.3226931  2.0035912 -0.6914797  1.2966556
##             Gfap    6.2252668 -0.8796777  1.6556304  0.6087384
##             Ednrb  10.9380415  4.5135260  3.0888286  1.6443596
##             
##              First four colnames:
##              240649020551054330404932383065726870513
##              274176126496863898679934791272921588227
##              323754550002953984063006506310071917306
##              87260224659312905497866017323180367450 
## 

6 Autocorrelation

Giotto Suite also now includes geospatial statistics and can perform global and local spatial autocorrelation

6.1 Global Statistics

viz = createSpatialNetwork(viz,
                           method = 'kNN',
                           name = 'kNN_network',
                           k = 8)

# returns as a data.table of results
res = spatialAutoCorGlobal(
  viz,
  method = 'moran',
  spatial_network_to_use = 'kNN_network',
  data_to_use = 'expression',
  expression_values = 'normalized'
)
## No spatial weight matrix found in selected spatial network
##  Generating distance matrix from kNN_network
print((res))
##       feat_ID         moran
##   1:     Mlc1  6.278703e-02
##   2:   Gprc5b  8.708937e-02
##   3:     Gfap  3.521745e-01
##   4:    Ednrb  1.488701e-01
##   5:     Sox9  1.229596e-01
##  ---                       
## 490: Blank-12 -2.874052e-03
## 491:    Ffar4 -1.280580e-03
## 492:    Qrfpr -3.748879e-03
## 493:     Chat -3.093600e-03
## 494: Blank-58 -7.066559e-05
print(res[moran == max(moran)])
##    feat_ID     moran
## 1: Slc17a7 0.5838803

6.2 Local Statistics

# local autocorrelation results are returned as a set of spatial enrichment data and saved in the giotto object
viz = spatialAutoCorLocal(viz,
                          method = 'moran',
                          data_to_use = 'expression',
                          expression_values = 'normalized')
## No spatial weight matrix found in selected spatial network
##  Generating distance matrix from kNN_network
## Attaching moran results as spatial enrichment: "expr_moran"
# local moran's I of the most autocorrelated gene according to the global results
spatPlot2D(viz,
           spat_enr_names = 'expr_moran',
           cell_color = res[moran == max(moran), feat_ID],
           color_as_factor = FALSE)

7 Generate Polygons

tx = getFeatureInfo(viz, feat_type = 'rna', return_giottoPoints = TRUE)
e = ext(tx) # get extent (spatial bounding box)

hexbin = tessellate(e, shape = 'hexagon', radius = 10, name = 'hexbin')
viz = setPolygonInfo(viz, hexbin)

pseudo_vis = makePseudoVisium(e)
viz = setPolygonInfo(viz, pseudo_vis, name = 'pseudo_visium')

viz = spatQueryGiottoPolygons(viz,
                              filters = list(pseudo_visium = 'all',
                                             hexbin = 'all'))
## Warning in spatQueryGiottoPolygons(viz, filters = list(pseudo_visium = "all", :
## NAs introduced by coercion

## Warning in spatQueryGiottoPolygons(viz, filters = list(pseudo_visium = "all", :
## NAs introduced by coercion
plot(hexbin)

plot(pseudo_vis)

plot(viz@spatial_info$query_poly)

8 GiottoData Package

GiottoData is a set of mini Giotto objects that are already put together and have some analyses run on them. These were created as spatial subsets of the full data making them more portable and easy to manipulate.

mini_viz = GiottoData::loadGiottoMini('vizgen')
mini_viz

9 Session Info

sessionInfo()
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] GiottoData_0.2.4  Giotto_3.3.1      data.table_1.14.8
## 
## loaded via a namespace (and not attached):
##  [1] reticulate_1.28     tidyselect_1.2.0    terra_1.7-39       
##  [4] xfun_0.39           bslib_0.4.2         listenv_0.9.0      
##  [7] lattice_0.20-45     colorspace_2.1-0    vctrs_0.6.2        
## [10] generics_0.1.3      htmltools_0.5.5     yaml_2.3.7         
## [13] utf8_1.2.3          rlang_1.1.1         R.oo_1.25.0        
## [16] jquerylib_0.1.4     pillar_1.9.0        glue_1.6.2         
## [19] withr_2.5.0         R.utils_2.12.2      lifecycle_1.0.3    
## [22] progressr_0.13.0    munsell_0.5.0       gtable_0.3.3       
## [25] R.methodsS3_1.8.2   future_1.32.0       codetools_0.2-18   
## [28] evaluate_0.21       labeling_0.4.2      knitr_1.42         
## [31] fastmap_1.1.1       parallel_4.2.1      fansi_1.0.4        
## [34] highr_0.10          Rcpp_1.0.11         backports_1.4.1    
## [37] checkmate_2.2.0     scales_1.2.1        cachem_1.0.8       
## [40] dbscan_1.1-11       jsonlite_1.8.4      parallelly_1.35.0  
## [43] farver_2.1.1        ggplot2_3.4.2       png_0.1-8          
## [46] digest_0.6.31       dplyr_1.1.2         grid_4.2.1         
## [49] rprojroot_2.0.3     scattermore_0.8     here_1.0.1         
## [52] cli_3.6.1           tools_4.2.1         magrittr_2.0.3     
## [55] sass_0.4.6          tibble_3.2.1        future.apply_1.10.0
## [58] pkgconfig_2.0.3     Matrix_1.5-4        rmarkdown_2.21     
## [61] rstudioapi_0.14     globals_0.16.2      R6_2.5.1           
## [64] igraph_1.4.2        compiler_4.2.1